HW 5.2

1 Helper Functions

Code

# https://github.com/wcmbishop/rayshader-demo/blob/master/R/image-size.R
define_image_size <- function(bbox, major_dim = 400) {
  # calculate aspect ration (width/height) from lat/long bounding box
  aspect_ratio <- abs((bbox$p1$long - bbox$p2$long) / (bbox$p1$lat - bbox$p2$lat))
  # define dimensions
  img_width <- ifelse(aspect_ratio > 1, major_dim, major_dim*aspect_ratio) %>% round()
  img_height <- ifelse(aspect_ratio < 1, major_dim, major_dim/aspect_ratio) %>% round()
  size_str <- paste(img_width, img_height, sep = ",")
  list(height = img_height, width = img_width, size = size_str)
}

# https://github.com/wcmbishop/rayshader-demo/blob/master/R/map-image-api.R
get_arcgis_map_image <- function(bbox, map_type = "World_Street_Map", file = NULL, 
                          width = 400, height = 400, sr_bbox = 4326) {
  require(httr)
  require(glue) 
  require(jsonlite)
  
  url <- parse_url("https://utility.arcgisonline.com/arcgis/rest/services/Utilities/PrintingTools/GPServer/Export%20Web%20Map%20Task/execute")
  
  # define JSON query parameter
  web_map_param <- list(
    baseMap = list(
      baseMapLayers = list(
        list(url = jsonlite::unbox(glue("https://services.arcgisonline.com/ArcGIS/rest/services/{map_type}/MapServer",
                                        map_type = map_type)))
      )
    ),
    exportOptions = list(
      outputSize = c(width, height)
    ),
    mapOptions = list(
      extent = list(
        spatialReference = list(wkid = jsonlite::unbox(sr_bbox)),
        xmax = jsonlite::unbox(max(bbox$p1$long, bbox$p2$long)),
        xmin = jsonlite::unbox(min(bbox$p1$long, bbox$p2$long)),
        ymax = jsonlite::unbox(max(bbox$p1$lat, bbox$p2$lat)),
        ymin = jsonlite::unbox(min(bbox$p1$lat, bbox$p2$lat))
      )
    )
  )
  
  res <- GET(
    url, 
    query = list(
      f = "json",
      Format = "PNG32",
      Layout_Template = "MAP_ONLY",
      Web_Map_as_JSON = jsonlite::toJSON(web_map_param))
  )
  
  if (status_code(res) == 200) {
    body <- content(res, type = "application/json")
    message(jsonlite::toJSON(body, auto_unbox = TRUE, pretty = TRUE))
    if (is.null(file)) 
      file <- tempfile("overlay_img", fileext = ".png")
    
    img_res <- GET(body$results[[1]]$value$url)
    img_bin <- content(img_res, "raw")
    writeBin(img_bin, file)
    message(paste("image saved to file:", file))
  } else {
    message(res)
  }
  invisible(file)
}

# https://github.com/wcmbishop/rayshader-demo/blob/master/R/elevation-api.R
get_usgs_elevation_data <- function(bbox, size = "400,400", file = NULL, 
                                    sr_bbox = 4326, sr_image = 4326) {
  require(httr)
  
  url <- parse_url("https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer/exportImage")
  res <- GET(
    url, 
    query = list(
      bbox = paste(bbox$p1$long, bbox$p1$lat, bbox$p2$long, bbox$p2$lat,
                   sep = ","),
      bboxSR = sr_bbox,
      imageSR = sr_image,
      size = size,
      format = "tiff",
      pixelType = "F32",
      noDataInterpretation = "esriNoDataMatchAny",
      interpolation = "+RSP_BilinearInterpolation",
      f = "json"
    )
  )
  
  if (status_code(res) == 200) {
    body <- content(res, type = "application/json")
    img_res <- GET(body$href)
    img_bin <- content(img_res, "raw")
    if (is.null(file)) 
      file <- tempfile("elev_matrix", fileext = ".tif")
    writeBin(img_bin, file)
    message(paste("image saved to file:", file))
  } else {
    warning(res)
  }
  invisible(file)
}

# https://github.com/wcmbishop/rayshader-demo/blob/master/R/find-image-coordinates.R
find_image_coordinates <- function(long, lat, bbox, image_width, image_height) {
  x_img <- round(image_width * (long - min(bbox$p1$long, bbox$p2$long)) / abs(bbox$p1$long - bbox$p2$long))
  y_img <- round(image_height * (1 - (lat - min(bbox$p1$lat, bbox$p2$lat)) / abs(bbox$p1$lat - bbox$p2$lat)))
  list(x = x_img, y = y_img)
}

2 Select Area of Interest: Honolulu, Hawaii

Code

# Bounding Box for Honolulu, Hawaii
bbox <- list(
  p1 = list(long = -158.30892000514868, lat = 21.723556651864115),
  p2 = list(long = -157.6240037655254, lat = 21.22700584076422)
)

# Confirming if bounding box is correct w/ leaflet
leaflet() %>%
  addTiles() %>% 
  addRectangles(
    lng1 = bbox$p1$long, lat1 = bbox$p1$lat,
    lng2 = bbox$p2$long, lat2 = bbox$p2$lat,
    fillColor = "transparent"
  ) %>%
  fitBounds(
    lng1 = bbox$p1$long, lat1 = bbox$p1$lat,
    lng2 = bbox$p2$long, lat2 = bbox$p2$lat,
  )

3 Define Image Size

Code

image_size <- define_image_size(bbox, major_dim = 600)

4 Get Image Overlay

Code

overlay_file <- "honolulu-map.png"
get_arcgis_map_image(bbox, map_type = "World_Imagery", file = overlay_file,
                     width = image_size$width, height = image_size$height, 
                     sr_bbox = 4326)
overlay_img <- png::readPNG(overlay_file)

5 Get Elevation Data

Code

# pull elevation data
elev_file <- file.path("elevation.tif")
get_usgs_elevation_data(bbox, size = image_size$size, file = elev_file,
                        sr_bbox = 4326, sr_image = 4326)

# load elevation data
elev_img <- raster::raster(elev_file)
elev_matrix <- matrix(
  raster::extract(elev_img, raster::extent(elev_img), buffer = 1000), 
  nrow = ncol(elev_img), ncol = nrow(elev_img)
)

zscale <- 30
# calculate rayshader water layer
watermap <- detect_water(elev_matrix)

6 Generate 3D Mappings

Code

n_frames <- 120 # increase n_frames from 100 to 120

# frame transition variables
min_depth_coeff <- 0.5
max_depth_coeff <- 2

waterdepthvalues <- mean(elev_matrix) * min_depth_coeff - (mean(elev_matrix) * (max_depth_coeff - min_depth_coeff)) * cos(seq(0, 2 * pi, length.out = n_frames)) ## adjust water depth range

thetavalues <- -90 + 45 * cos(seq(0, 2*pi, length.out = n_frames))

# labels
label1 <- list(text = "Honolulu")
label1$lon <- -157.8581
label1$lat <- 21.3099

label2 <- list(text = "Pearl City")
label2$lon <- -157.9715
label2$lat <- 21.3910

# generate .png frame images
img_frames <- paste0("drain", seq_len(n_frames), ".png")
for (i in seq_len(n_frames)) {
  message(paste(" - image", i, "of", n_frames))
  elev_matrix %>%
    sphere_shade(texture = "imhof4") %>%
    add_water(watermap, color = "imhof3") %>%
    add_overlay(overlay_img, alphalayer = 0.5) %>%
    plot_3d(elev_matrix, solid = TRUE, shadow = TRUE, zscale = zscale, 
            water = TRUE, watercolor = "imhof3", wateralpha = 0.8, 
            waterlinecolor = "#ffffff", waterlinealpha = 0.5,
            waterdepth = waterdepthvalues[i]/zscale, 
            theta = thetavalues[i], phi = 45, zoom = 0.65,
            fov = 60)
  
  # find relative image coordinates based on lon, lat of labels
  label1$pos <- find_image_coordinates(label1$lon, label1$lat, bbox, image_size$width, image_size$height) 
  label2$pos <- find_image_coordinates(label2$lon, label2$lat, bbox, image_size$width, image_size$height)  

  render_label(elev_matrix, x = label1$pos$x, y = label1$pos$y, z = 500, 
             zscale = zscale, text = label1$text, textsize = 1, linewidth = 3)
  render_label(elev_matrix, x = label2$pos$x, y = label2$pos$y, z = 500, 
             zscale = zscale, text = label2$text, textsize = 1, linewidth = 3)
  render_snapshot(img_frames[i])
  rgl::clear3d()
}

7 Convert 3D Mappings into GIF

Code

# build gif
magick::image_write_gif(magick::image_read(img_frames), 
                        path = "water_rise.gif", 
                        delay = 12/n_frames) # increase delay for smoother animation